function F = RF_beamformer_design(Mx, My, d, AoD_groups)
% RF_BEAMFORMER_DESIGN RF波束成形设计（Algorithm 2）
%   输入：
%     Mx, My: URA尺寸（x/y方向天线数）
%     d: 天线间距（波长倍数）
%     AoD_groups: [1×G]结构体数组（extract_AoD_support输出）
%   输出：
%     F: [M×b]RF波束成形矩阵，M=Mx*My，b=总波束数
%
%   文献依据：公式(19)-(20)

    % ========== 1. 参数验证 ==========
    validate_parameters(Mx, My, d, AoD_groups);
    
    % ========== 2. 初始化输出矩阵 ==========
    M = Mx * My; % 总天线数
    F = []; % 初始化RF矩阵
    
    % ========== 3. 生成量化角度网格 ==========
    lambda_x = linspace(-1, 1, Mx); % x轴量化角度
    lambda_y = linspace(-1, 1, My); % y轴量化角度
    
    % ========== 4. 遍历每个用户组 ==========
    for g = 1:length(AoD_groups)
        % 4.1 获取当前组AoD支持区域
        AoD = AoD_groups(g);
        
        % 4.2 筛选覆盖AoD区域的量化角度对（公式19）
        selected_pairs = [];
        for i = 1:Mx
            for j = 1:My
                % 计算当前量化角度对的边界
                lambda_x_min = lambda_x(i) - 1/Mx;
                lambda_x_max = lambda_x(i) + 1/Mx;
                lambda_y_min = lambda_y(j) - 1/My;
                lambda_y_max = lambda_y(j) + 1/My;
                
                % 检查是否覆盖AoD区域（公式19）
                if is_covering_AoD(lambda_x_min, lambda_x_max, ...
                                  lambda_y_min, lambda_y_max, AoD)
                    selected_pairs = [selected_pairs; i, j];
                end
            end
        end
        
        % 4.3 构建当前组的RF波束成形矩阵（公式20）
        F_g = zeros(M, size(selected_pairs, 1));
        for idx = 1:size(selected_pairs, 1)
            i = selected_pairs(idx, 1);
            j = selected_pairs(idx, 2);
            
            % 生成导向矢量
            e_vec = steering_vector(lambda_x(i), lambda_y(j), Mx, My, d);
            F_g(:, idx) = e_vec;
        end
        
        % 4.4 添加到总矩阵
        F = [F, F_g];
    end
    
    % ========== 5. 正交性验证 ==========
    validate_orthogonality(F);
end

%% ========== 子函数：导向矢量生成 ==========
function e_vec = steering_vector(gamma_x, gamma_y, Mx, My, d)
% 生成导向矢量（公式3）
    % x方向导向矢量
    phase_x = -2*pi*d*gamma_x*(0:Mx-1)';
    e_x = exp(1i*phase_x) / sqrt(Mx);
    
    % y方向导向矢量
    phase_y = -2*pi*d*gamma_y*(0:My-1)';
    e_y = exp(1i*phase_y) / sqrt(My);
    
    % Kronecker积
    e_vec = kron(e_x, e_y);
end

%% ========== 子函数：AoD覆盖检查 ==========
function flag = is_covering_AoD(lambda_x_min, lambda_x_max, ...
                                lambda_y_min, lambda_y_max, AoD)
% 检查量化角度对是否覆盖AoD区域（公式19）
    % 提取AoD边界
    gamma_x_min = AoD.gamma_x_range(1);
    gamma_x_max = AoD.gamma_x_range(2);
    gamma_y_min = AoD.gamma_y_range(1);
    gamma_y_max = AoD.gamma_y_range(2);
    
    % 检查覆盖条件
    flag_x = (gamma_x_min >= lambda_x_min) && (gamma_x_max <= lambda_x_max);
    flag_y = (gamma_y_min >= lambda_y_min) && (gamma_y_max <= lambda_y_max);
    
    flag = flag_x && flag_y;
end

%% ========== 子函数：参数验证 ==========
function validate_parameters(Mx, My, d, AoD_groups)
    % 天线参数
    if ~isscalar(Mx) || Mx < 1 || mod(Mx,1) ~= 0
        error('Mx must be a positive integer');
    end
    if ~isscalar(My) || My < 1 || mod(My,1) ~= 0
        error('My must be a positive integer');
    end
    if ~isscalar(d) || d <= 0
        error('d must be a positive scalar');
    end
    
    % AoD_groups验证
    if ~isstruct(AoD_groups) || isempty(AoD_groups)
        error('AoD_groups must be a non-empty structure array');
    end
    required_fields = {'gamma_x_range', 'gamma_y_range'};
    for i = 1:length(required_fields)
        if ~isfield(AoD_groups, required_fields{i})
            error('AoD_groups missing required field: %s', required_fields{i});
        end
    end
end

%% ========== 子函数：正交性验证 ==========
function validate_orthogonality(F)
% 验证F^H F ≈ I（文献公式15）
    if isempty(F)
        return;
    end
    
    % 计算正交误差
    orth_error = norm(F'*F - eye(size(F,2)), 'fro');
    
    % 设置容差阈值
    tol = 1e-6;
    
    if orth_error > tol
        warning('RF波束成形矩阵正交性误差: %.2e', orth_error);
    end
end